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ABSTRACT 



Aims. We exploit the deepest existing far-infrared (FIR) data obtained so far by Herschel at 100 and 160 pm in the GOODS-N, as part 
of the PACS Evolutionary Probe (PEP) survey, to derive for the first time the evolution of the rest-frame 60-pm, 90-/jm, and total IR 
luminosity functions (LFs) of galaxies and AGNs from z = to unprecedented high redshifts (z~2—3). 

Methods. The PEP LFs were computed using the \fW max method. The FIR sources were classified by means of a detailed broad- 
band SED-fitting analysis and spectral characterisation. Based on the best-fit model results, k-correction and total IR (8-1000 pm) 
luminosity were obtained for each source. LFs (monochromatic and total) were then derived for various IR populations separately in 
different redshift bins and compared to backward evolution model predictions. 

Results. We detect strong evolution in the LF to at least z~2. Objects with SEDs similar to local spiral galaxies are the major 
contributors to the star formation density (SFD) at zl03, then, as redshift increases, moderate SF galaxies - most likely containing 
a low-luminosity AGN - start dominating up to z-1.5. At z ~ 1.5 the SFD is dominated by the contributions of starburst galaxies. In 
agreement with previous findings, the comoving IR LD derived from our data evolves approximately as (1 + z) 3 " 8±0 " 3 up to z~l, there 
being some evidence of flattening up to z~2. 

Key words. Galaxies: evolution - Galaxies: active - Galaxies: starburst - Cosmology: observations - Infrared: galaxies 



1. Introduction 

Extragalactic surveys in the infrared (IR) represent a key tool 
for understanding galaxy formation and evolution, the luminos- 
ity function (LF) computed at different redshifts being the most 
direct method for exploring the evolution of a galaxy population. 
Deep surveys with IRAS and ISO allowed the first studies of the 
IR-galaxy LF at z! 0.3 (Saunders et al. 1990) and atz! 1, respec- 
tively (Pozzi et al. 2004). With Spitzer 24-yum data, it was pos- 
sible to study the evolution of the mid-infrared (MIR) LF up to 
z-2.5 (e.g. Le Floc'h et al. 2005; Caputi et al. 2007; Rodighiero 
et al. 2010), while, even with the deepest Spitzer 70-//m data, 
only z~ 1 could be reached in the far-infrared (FIR) (Magnelli et 
al. 2009). Ground-based observations in the mm/sub-mm range, 
probing the evolution of the most distant (z ! 2) and luminous 
dusty galaxies, have so far been limited to the identification of 
sources at the very bright end of the luminosity function (e.g., 
Chapman et al. 2005). All of these works detected strong evolu- 
tion in both luminosity and density, indicating that IR galaxies 
were more numerous and more luminous in the past. The de- 
rived rates of evolution of the IR sources are significantly higher 
than observed at any other wavelength. Since the rest-frame IR 
spectral energy distributions (SED) of star-forming galaxies and 
AGNs peak at 60-200 //m, to measure their bolometric lumi- 
nosity and evolution with z we need to observe in the FIR/sub- 



mm regime. However, the detection of large numbers of high-z 
sources at the peak of their IR SED could not yet been achiev- 
able, due to source confusion and/or low detector sensitivity. 

In the present paper, we discuss observations performed 
with the PACS instrument (Poglitsch et al. 2010) onboard the 
Herschel Space Observatory (Pilbratt et al. 2010). Herschel is 
the first telescope allowing us to detect the FIR population to 
high redshifts (z~3-4) and to derive its rate of evolution by a de- 
tailed LF analysis. The PACS Evolutionary Probe (PEP, PI D. 
Lutz) Surveys are designed specifically to this purpose. With the 
observations dedicated to the GOODS-N by the PEP Science 
Demonstration Program (SDP), we first address the main goal 
of determining the evolution with redshift of the galaxy and 
AGN LF in the FIR domain. The availability of a deep multi- 
wavelength catalogue in the GOODS-N is crucial for analysing 
the SEDs, obtaining k-corrections and total IR luminosities, and 
studying the LF and its evolution separately for the different IR 
populations. The measure of the total IR luminosity allows us to 
derive for the first time the total IR LF and its evolution directly 
from FIR data (see e.g. Nordon et al. (2010), this issue, for de- 
scriptions of the failure, at least at z>1.5, of previous total IR 
luminosity extrapolations from the MIR). 

Throughout this paper, we assume a ACDM cosmology with 
H = 71 km s- 1 Mpc- 1 , Q m = 0.27, and Q. A = 0.73. 



* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im- 
portant participation from NASA. 
** e-mail: carlotta.gruppioni@oabo.inaf.it 



2. Multiwavelength source characterisation 

Our reference samples were the PEP SDP blind catalogues at 
100 and 160 pm to the 3<x level (corresponding to ~3 and 5.7 
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mJy, respectively) in the GOODS-N. The PEP sources were 
cross-correlated with the 24-/mi and multi-wavelength cata- 
logues by means of a likelihood ratio technique (Sutherland & 
Saunders 1992). We refer to the Appendix of Berta et al. (2010) 
for a detailed description of the data and the source identifica- 
tion. To maximise the fraction of identifications, we limited our 
study to the area covered by the ACS (~150 arcmin 2 ), obtaining 
227 and 253 sources at 100 and 160 /mi, respectively (216 and 
237, respectively, with either spec- (~70%) or photo-z). 

We used of all the available data to derive the SEDs of our 
PEP sources, which we interpreted and classified by perform- 
ing a fit (using the Le PHARE software; Ilbert et al. 2006) with 
a semi-empirical template library (Polletta et al. 2007) repre- 
sentative of different classes of IR galaxies and AGNs (spirals, 
starbursts (SB), QSOs, Seyfert2/1.8, and composite AGN+SB 
objects). Three main SED classes are found to reproduce most 
of our sources: spiral, SB, and AGN2 (Seyfert2/1.8). The spiral 
SEDs show no clear signs of enhanced SF or nuclear activity 
(see the online Fig. [5}, the FIR bump being characterised by rel- 
atively cold dust (T</„„~20 K). On the other hand, SB templates 
are characterised by warmer (T</ M;tf ~40^-5 K), more pronounced 
FIR bumps and significant UV extinction, indicative of intense 
SF activity. Templates of moderately luminous AGNs (repre- 
senting Seyfert2/1.8 galaxies) are characterised by a "flattening" 
in the 3-10 /mi spectrum (suggesting detection of an AGN in the 
wavelength range where the host galaxy SED has a minimum) 
and a FIR bump dominated by SF and intermediate (in terms of 
both energy and T^ 1(J ,) between spirals and SBs. Although they 
can be considered as SF galaxies at the wavelengths relevant 
to this work, we prefer to refer to them as "AGN2" throughout 
the paper, to keep in mind that they probably contain an AGN, 
whose presence, though not dominant in the FIR, might be very 
important for analysis in other bands (i.e., in the X-rays or the 
MIR). In general, the considered templates provide very good 
fits to the SEDs of our PEP sources. However, in ~10— 12% of 
cases the observed SEDs are very well reproduced by the tem- 
plates over the entire UV/optical/NIR/IRAC range, while they 
are systematically underestimated in the MIPS/PACS range. In 
these cases, the PEP 100- and 160-/im data in flux density are 
always higher by up to a factor of ~4 than the template at the 
same A's. This happens mainly for the Seyfert2/1.8 templates 
(for about 40% of the PEP sources fitted by these SEDs) and in 
less frequent cases also for the spiral ones. We therefore con- 
structed three new templates with a rest-frame 0.1-15 /mi spec- 
trum similar to that of Polletta et al. (2007), but with a higher 
FIR bump, by averaging together (in /l-bins) the observed rest- 
frame SEDs (normalised to K s band) exhibiting an excess in the 
FIR and fitted by the same template (Seyfert2/1.8/Sdm; see an 
example in Fig.|6]of the online edition). We also added to the li- 
brary three new SB templates obtained by interpolating between 
the Sdm and the SB SED NGC6090 (to fill the large gap ex- 
isting in the library between spirals and SBs). The final results 
of the SED-fitting are reported in the online Table [T] Our SED- 
fitting results show a clear trend with z, with the sources becom- 
ing "warmer", "more active", and "more IR luminous" with in- 
creasing z (see the online Fig. 17] where the 160-//m/100-//m flux 
ratio versus z is plotted for the PEP 160-/mi selected sources). 
In particular, spiral galaxies are found exclusively at low-z (<1), 
dominating at z<0.4 and coexisting with AGN2 and SB galax- 
ies at 0.4<z<l. The AGN2 are observed up to z~ 1.5-2, with a 
strong concentration at 0.4<z<l. The SB galaxies are found at 
any redshift in the range z~0.4-3.0 and dominate between z^l 
and z-2.5. The composite objects and the QSOs are found ex- 
clusively at high-z (>1.5 and >3, respectively). Although sample 
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Fig. 1. Rest-frame LF at 60 /mi (top) and 90 /mi (bottom) in dif- 
ferent redshift bins, computed as \fW max (filled circles). The grey 
dashed line at 60 /mi is the Saunders et al. (1990) LLF, while 
the diagonal crosses at 90 /mi represent the LLF of Serjeant et 
al. (2004). The GP2010 evolutionary model predictions are also 
shown as a black solid line. 



selection effects may also play a role in the observed trend with 
z, we conclude that sources with enhanced SF activity dominate 
at higher-z, gradually becoming "colder","less active", and "less 
IR luminous" as we move towards lower redshift. Similar re- 
sults were found by Gruppioni et al. (2008), who performed the 
same SED-fitting study presented here, but for a MIR-selected 
sample. The results of Gruppioni et al. (2008) were then used 
by Gruppioni & Pozzi (2010, in preparation, hereafter GP2010) 
to construct a backward evolution model of IR sources, assum- 
ing separate evolutionary paths for the five populations described 
above. 



3. Luminosity functions 

The sizes of the PEP samples are sufficiently large to allow a di- 
rect determination of the FIR LF to z~2, and provide constraints 
on evolution to z~3-4. For the median redshift of the 100- and 
160-//m PEP samples (~0.68 and 0.82 respectively), we com- 
puted the rest-frame LFs at 60 and 90 //m, using the 1/V max 
method (Schmidt 1968). We divided the two samples into six 
redshift bins (slightly different for the two wavelengths) selected 
to be almost equally populated (~50 sources) at least up to z~2. 
In the two higher, less populated, redshift bins, the characteris- 
tic luminosity L* is not constrained by data, but some hints of 
evolution can be detected (note that at z>3 most of the objects 
do have power-law SEDs and only photo-z's, which are very un- 
certain due to the SED shape). The completeness and spurious 
correction factors given by Berta et al. (2010) were applied to 
our PEP data. The redshift incompleteness does not affect our 
conclusions, since ~95% of our sources have a redshift. The re- 
sults of the computation of our 60- and 90-/mi LFs are shown in 
Fig.[T]and reported in Tables|2]and[3]of the online edition of this 
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Fig. 2. Rest-frame 90-yum LF at different redshifts, as com- 
puted from PEP (coloured filled triangles) and predicted 
by the GP2010 model for the different IR populations 
(green dashed line: spirals (three top panels); red dot- 
dot-dot-dashed: AGN2 (second row); cyan dot-dashed: SBs 
(third row and bottom left panel); magenta dotted: composite 
(bottom middle); blue long-dashed: AGN1 (bottom right)). The 
total observed (open circles) and predicted (continuous line) LFs 
of Fig.[T](Z?o£ torn) are reported. The vertical dashed lines are the 
luminosities above which we expect our survey to be complete 
in that z-bin for that population. 
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Fig. 3. Total IR LF from our PEP data (black filled circles). 
Other results from the literature are plotted for comparison (di- 
agonal crosses: LLF of Sanders et al. (2003); red filled trian- 
gles: Vaccari et al. (2010); green filled triangles: Le Floc'h et 
al. (2005); blue filled squares: Rodighiero et al. (2010); orange 
filled squares: Chapman et al. (2005)). The black solid, green 
dot-dashed, blue dotted, and red long-dashed lines represent the 
expectations from the GP2010, Lagache et al. (2004), Valiante 
et al. (2009), and Le Borgne et al. (2009) models, respectively. 



paper. The vertical dashed lines in each z-bin are the luminosi- 
ties below which we expect our samples to be incomplete, given 
that at fainter luminosities not all galaxy types are observable 
(depending on their SED; Ilbert et al. 2004). For comparison, 
we overplot the LLFs at 60 pm and 90 pm from S aunders ( 1 990) 
and Serjeant et al. (2004), respectively, and the predictions of 
the GP2010 backward evolution model. Strong luminosity evo- 
lution is necessary to explain the observed LF, even in the lower 
z-bin (z<0.4-0.5). To investigate the different evolutionary paths 
of the various IR populations, we computed the 1/V mfl x LF sepa- 
rately for the five classes defined by the SED-fitting analysis. In 
Fig. [2] we show the 90-pm LFs derived from the 160-//m sample 
for the different SED classes, plotted in the most representative 
z-bins for each of them. For comparison, we also show the sin- 
gle population LFs of the GP2010 model. The close agreement 
between data and model implies that both the evolution and the 
template SEDs (modified in the FIR as described in Sect. |2]i as- 
sumed by the model (based on MIR results at z! 1.5) to represent 
the different IR classes can also be considered a fair representa- 
tion of the FIR populations, at least up to z~3-4. In agreement 
with the change in population with z implied by the SED-fitting 
analysis and discussed in Sect. |2l we observe that at lower red- 
shifts (z<0.5) the LF is dominated by galaxies with spiral SEDs, 
while at 0.5<z<0.8 Seyfert2/1.8-SED sources appear to be more 
important than the normal galaxy population, becoming domi- 
nant around L* (10 105 L G <L9o<10 n 5 L ) at 0.5<z<1.2. The SB 
population prevails at the higher luminosities (L9o>10 1L5 L G ) at 
any redshifts, dominating also the L* regime at zl2. The com- 
posite and quasar populations are never dominant in the FIR, 
contributing to the 60- and 90-jt/m LFs mainly at z>2. 



3. 1 . The total IR luminosity function 

We integrate the best-fit SED of each source at 8</l rejf <1000 
pm to derive the total IR luminosities (L/«=L[8-1000 pm]). Our 
approach is similar to that of other studies based on MIR se- 
lected galaxy samples (i.e. Le Floc'h et al. 2005; Rodighiero et 
al. 2010), but this is the first time that the SEDs are accurately 
constrained by sufficiently deep data in the FIR domain and not 
simply extrapolated from the MIR to the FIR/sub-mm. The total 
IR LF up to z~3 was computed (values reported in the online 
Table HI and compared with other derivations available in the 
literature, as shown in Fig. [3] The total IR LF of Sanders et al. 
(2003) is plotted as a local reference, in addition to the predic- 
tions of four backward evolution models as references of evolu- 
tion. Globally, data and model expectations seem to agree well 
over the whole z-range. The agreement between the PEP total 
IR LF and those computed from 24-jum (ECDFS 0.3<z<1.2: Le 
Floc'h et al. 2005; VVDS+GOODS 0<z<2.5: Rodighiero et al. 
2010) and SPIRE surveys (Hermes SDP: 0<z<0.2; Vaccari et al. 
2010) is very good in the z-bins where the latter are available. 
The total IR LF of sub-mm galaxies from Chapman et al. (2005) 
at z~2.5 represents very well the bright end of the total IR LF, 
but seems to be affected by significant incompleteness at z~0.9. 
The sub-mm surveys clearly provide the most reliable probe of 
the evolution of the more distant and more luminous dusty galax- 
ies, but are not sensitive enough to sample the low-intermediate 
redshift Universe. Our PEP total IR LF evolves in luminosity as 
~(1 + z) 4 ' 1 - 02 up to z=1.5±0.3 (though considerable degeneracy 
is found between strict luminosity evolution and a combination 
of both density and luminosity). At 1.5 tzl 2.5-3, the evolution 
rate appears to remain almost constant. 
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Fig. 4. Evolution of the total IR LD (or SFD) with z. The results 
of integrating our observed total IR LF in each z-bin are shown as 
black filled circles. The 3cr best-fitting envelope to the Hopkins 
& Beacom (2006) data collection from different surveys (yel- 
low area), the results of Le Floc'h et al. (2005) and Rodighiero 
et al. (2010) as total IR LD up to z~l (green area) and z~2.5 
(blue filled squares), respectively, and the contribution to the 
SFD of radio detected sub-mm sources (orange filled squares) 
from Chapman et al. (2005) are also reported. The differently 
coloured lines are the GP2010 model predictions in terms of total 
IR LD (black solid), IR LD to the PEP luminosity completeness 
limit (grey solid), and single contributions from the different IR 
populations (same as in Fig.|2]i. 



3.2. The comoving IR luminosity density 

One of the main reasons for computing the total IR LF is that 
it allows a direct estimate to be made of the total comoving IR 
luminosity density (LD) and the star formation density (SFD) as 
a function of z, which are crucial tools for understanding galaxy 
formation and evolution. We convert L/« into star formation rate 
(SFR) by using the relation (for a Salpeter IMF) SFR [M Q yr _I ] 
= 1.7X10" 10 L IR [Lq] (Kennicutt, 1998), then compute the SFD 
directly from our observed PEP LF, as well as from the GP2010 
model. Since the model provides a fair representation of the data, 
the predicted SFD can be interpreted as what we could expect 
from an "ideal" FIR survey (or collection of FIR surveys) cov- 
ering the whole range of luminosities in each redshift bin. Our 
data are able to give a robust estimate of the LD only at z 1 1-1 .2, 
providing just lower limits at higher-z, given that with increas- 
ing redshift our total IR LF starts to become incomplete at the 
lower luminosities, the completeness limit moving towards the 
bright end of the luminosity distribution (at z>l, we sample just 
L/^>L* S ). The IR LD (or SFD) estimates from either data or 
model are shown in Fig. [4] To model the luminosity incomplete- 
ness of our data, we integrated the GP2010 total IR LF in each 
z-bin down to the completeness limit shown in Fig. [3] obtaining 
a SFD prediction (grey solid line of Fig. HI that is totally con- 
sistent with our data estimates. For comparison, we have also 
shown the SFD obtained from different surveys by Hopkins & 
Beacom (2006) and the results of 24-/im surveys by Le Floc'h 
et al. (2005) and Rodighiero et al. (2010) and sub-mm surveys 
by Chapman et al. (2005). We find very close agreement with 
both previous results, based on either IR data or data at differ- 
ent wavelengths, and the GP2010 model predictions up to z~l. 
From PEP data, we find that the comoving IR LD evolves as 
~(1 + z) 3 ,8±0 ' 3 up to z~l, The model predicts a predominance of 
normal spiral galaxies only at very low redshifts (z!0.3), when 
sources with Seyfert2/1.8 SEDs begin to dominate the SFD up to 



z~1.5. The SB galaxies then become the prevalent population up 
to the highest redshifts. We note that we have not subtracted the 
AGN populations from our SFD calculation, since the FIR part 
of the spectrum of all the templates (but the very few AGN1) 
is probably dominated by SF activity. Even in the Seyfert2/1.8 
templates, the AGN activity dominates over the stellar emission 
only at A res ,~3-IQ fim. Given the power-law shape of their spec- 
tra, it is less clear at which A the AGN should dominate the light 
of the composite SEDs, and consequently what fraction of the to- 
tal IR luminosity is related to accretion activity rather than SF in 
these objects. However, these populations are never predominant 
in our FIR survey and even a small fraction of contamination re- 
lated to accretion activity occurring in these objects (mainly at 
high-z) will not significantly affect our results. 

4. Conclusions 

We have used of the deepest Herschel 100- and 160-yum data 
from the PEP SDP survey in the GOODS -N to characterise the 
evolution of the galaxy and AGN FIR LF and LD across the red- 
shift range Qtzt 3. To summarise, in the present work we have: 

1 . completely characterised the multiwavelength SEDs of the 
PEP sources by performing a detailed SED-fitting analysis, 
tightly constraining (and in some cases opportunely modify- 
ing) the templates at long wavelengths. 

2. computed the first rest-frame LFs at 60 and 90 /urn up to z~4, 
finding that a significant amount of evolution is required to 
reproduce the LF at least up to z~2, the different SED-classes 
exhibiting different evolutionary behaviours. 

3. integrated the SEDs at A rest =8-1000 fim and computed the 
total IR LF up to z~3, finding luminosity evolution oc(l + 
z) -o.2 up to z=1.5 ± 0.3, then an almost constant evolution 
rate toz~2.5-3. 

4. derived the evolution of the comoving total IR LD (and 
SFD), which was found to increase as (1 + z ) 3 - 8±a3 up to 
Z~l, Spiral galaxies are the main contributors to the SFD at 
z;0.3, moderate SF galaxies (which most likely harbour a 
low luminosity AGN) becoming increasingly important up 
to z-1.5. At this redshift, the SFD is representative of SB 
galaxies. 
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Table 1. SED-fitting results 



AQxm) 


Spiral 


Starburst 


Composite 


AGN2 


AGN1 


100 
160 


60 (28%) 
63 (27%) 


70 (32%) 
82 (35%) 


13 (6%) 
13 (5%) 


71 (33%) 
76 (32%) 


2 (1%) 

3 (1%) 



Table 2. Rest-frame 60 fim luminosity function in the GOODS-N 



log 


C>L(,0imt (Lq) 






log 10 O (Mpc 


-' dex" 1 ) 










0.0<z<0.4 


0.4<z<0.6 


0.6<z<l.l 


l.l<z<1.9 


1.9<z<2.9 


2.9<z<4.2 




8.3 


-2.03 ± 0.73 














8.9 


-1.89 ±0.41 














9.5 


-2.60 ± 0.34 














10.1 


-2.79 ±0.25 


-2.72 ±0.39 












10.7 


-3.47 ± 0.45 


-2.73 ±0.17 


-2.97 ± 0.50 










11.3 




-3.60 ±0.35 


-3.16 + 0.14 - 


-3.84 ±0.43 








11.9 






-4.18 ±0.33 - 


-3.96 ± 0.26 


-3.98 ±0.81 






12.5 








-5.06 ±0.58 


-4.54 ±0.38 






13.1 










-5.68 ±0.89 


-5.39 ±0.72 



Table 3. Rest-frame 90 fim luminosity function in the GOODS-N 



log 


oLgo/mt (Le) 






logioO (Mpc - ^ dex -1 ) 










0.0<z<0.5 


0.5<z<0.8 


0.8<z<1.2 1.2<z<2.0 


2.0<z<3.0 


3.0<z<4.2 




9.2 


-1.63 ±0.41 












9.8 


-2.44 ± 0.33 












10.4 


-2.64 ±0.21 


-2.28 ± 0.40 










11.0 


-3.20 ±0.25 


-2.95 ±0.21 


-2.86 ± 0.26 








11.6 




-3.92 ±0.38 


-3.45 ±0.19 -3.18 ±0.47 








12.2 






-4.80 ±0.71 -4.37 ±0.29 


-4.27 ± 0.35 


-4.31 ±0.61 




12.8 








-5.19 ±0.58 


-5.29 ±0.71 
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Fig. 5. Template SEDs of a Starburst (M82: cyan dot-dashed line), a Seyfert 2 (red dot-dot-dot-dashed) and a Sdm galaxy (green 
dashed) from the Polletta et al. (2007) library, normalised to the fi-band flux. 
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Fig. 6. Rest-frame observed SEDs of the PEP 100- and 160-/im sources fitted by the Polletta et al. (2007) Seyfert 2 template (orange 
dashed line) with underestimated FIR. The red solid line shows the new "FIR Seyfert 2" template obtained by averaging the observed 
data in wavelength bins. 
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Fig. 7. 160-pm/lOO-jt/m flux ratio versus z both observed for the differently classified PEP 160-/im sources (symbols) and derived 
from the characteristic templates (curves): green circles and dashed line represent the normal galaxy population; cyan stars and dot- 
dashed line are for SB; red triangles and dot-dot-dot-dashed line are for AGN2; magenta squares and dotted line are for composite 
objects; blue triangles and long-dashed line are for AGN1. The upward arrows are for sources not detected at 100 fim (for which 
the 3 mJy flux limit has been considered). 



